Sn63 


*■ 

» * 

3&P • 


V 


IMPROVED TWO -DIMENSIONAL -KINETICS 


COMPUTER PROGRAM 


. NAS8-35931 


1 


Third Quarterly P 
July 1984 through 


ingress Repor 
30 September 


t 

1984 


/ 


A/fe# -c/r- /V//zy 

(NASA-CR-4 7Q7 -fr») IMPROVED 

TWO- DIMENSION AL-K3 NETICS COMPUTER PROGRAM 
Quarterly Progress Report, 1 Jul - 30 Sep. 

1S64 (Software and Engineering Associates, 

Inc.) 23 p CSCL 09B G3/61 


N87-1 1508 


Unclas 

43852 


Prepared by: Software and Engineering Associates, Inc. 

1560 Brookhollow Drive, Suite 203 
Santa Ana, California 92705 


1 



Prepared 


Prepared 


Sn63 


IMPROVED TWO-DIMENSIONAL-KINETICS 
COMPUTER PROGRAM 


NAS8-35931 


Third Quarterly Progress Report 
1 July 1984 through 30 September 1984 


for: George C. Marshall Space Flight Center 

Marshall Space Flight Center 
Alabama 35812 


by: Software and Engineering Associates, Inc. 

1560 Brookhollow Drive, Suite 203 
Santa Ana, California 92705 


1 


IMPROVED TWO-DIMENSIONAL-KINETICS COMPUTER PROGRAM 


1 . 0 BACKGROUND 


Future Orbit Transfer Vehicles (OTV) presently under 
consideration need rocket engines delivering a high specific 
impulse. This high performance can be obtained with large 
area ratio thrust chambers using oxygen with either hydrogen 
or hydrocarbon fuels. In the projected nozzles the 
combustion products are expanded to low pressure and 
temperature levels at high Mach-number flow, a domain which 
has not been experienced with existing rocket engines. 

Some modifications have recently been incorporated in 
the T wo-D imens i onal- Ki net i cs (TDK) computer program, such 
as: Condensed phase flow simulation, treatment of shock 
waves induced by the wall curvature, and the coupling of the 
TDK code with the Boundary Layer Module (BLM) for an 
improved thrust chamber specific impulse prediction. 

Recent boundary layer calculations for large area ratio 
nozzles have revealed that the boundary layer becomes very 
thick in the region of high Mach-numbers . This result is 
quite different from rocket engines built up to now and 
requires an examination of the existing thrust ' loss 
calculation method due to the viscous effects adjacent to 
the wall. Furthermore, the knowledge of the Knudsen-number 
is desirable which differentiates between continuum, slip 
and free molecular flow, and thereby identified when the 
analysis, presently based on Newtonian fluid flow, becomes 
questionable or invalid. 
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In order to obtain highly accurate performance 
predictions, the best supporting data for equilibrium and 
finite rate chemistry must be prepared. Furthermore, the 
large area ratio nozzles require many more calculation steps 
to cover the complete nozzle flow field. This increases 
computation time significantly and favors an error 
accumulation which should be reduced as much as possible. 
Another point of importance connected with the OTV is the 
existence of ionized combustion products which permits 
vehicle tracing or interferes with the transfer of 
communication or command signals. 

The current effort of shock wave modeling induced by 
the curvature of the wall should be expanded to permit the 
shock simulation caused by wall contour discontinuities. 
Also the Mach shock existence in the center of the flow 
field interacting with other shocks as well as the 
analytical treatment of multiple shock waves may 
significantly affect the interior nozzle flow and associated 
performance and need to be simulated. 
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2.0 OBJECTIVE 


The objective of this effort is to increase the 
analytical capability of the existing TDK/BLM computer 
program for performance of OTV thrust chambers. Areas which 
need further examination and improvement are the thick 
boundary layer inviscid core flow interaction and the 
related thrust loss calculation. To warrant highly accurate 
results the provision of the best available program input 
data is mandatory as well as the use of sophisticated 
modeling techniques to reduce computation time with advanced 
error control criteria. The simulation of wall shocks, 
Macn-shocks in addition to shocks induced by iarge concave 
wall curvature is essential since this interaction with each 
other produces flow field changes which in turn affect the 
nozzle performance. 
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3 . 0 WORK PERFORMED DURING THE REPORTING PERIOD 

The program plan for this work is presented in Figure A. 
Progress in each of the work tasks during the reporting period 
is presented in the sub-paragraphs 3.1 through 3.5. 


3 . 1 TDK-BLM INTERFACE 


Interaction between the inviscid core flow and the wall 
boundary layer is expressed through 1st order by the definition 
of displacement thickness, 6*. We have rederived <$* and 9 for 
the nozzle wall boundary layer situation for the purpose of 
obtaining a more accurate representation, especially when the 
boundary layer is thick. In particular the radial coordinate, 
r, has not been factored. Integral expressions are obtained 
that are somewhat different than the classical expressions. For 
example, the expression fo^ 5* is found to be: 
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(1 
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r +r 
1 e 


pur, 
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This is an implicit expression, since conditions at e are 
displaced from the wall by a distance' 6*. The subscript 1 
refers to a displacement equal to the velocity thickness. It is 
found that the boundary layer thrust deficit does not use this 
integral, but a somewhat different integral, as follows: 


* f / \ — p+r 

<s h - (.1 - pur J dy , r - 1 w . 

b J 0 1 w “2“ 

p e u e r 1 w 

This is a explicit expression. The subscript w refers to 
conditions at the wall. 

Similar remarks apply to the momentum thickness, 0. 


FIGURE PROGRAM PLAN (as of 30 Sept 1984) 
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3.2 CHEMISTRY IMPROVEMENTS 


GRKif: v.;;.; • ^ 

OF POOR QUALITY 


Task 3.2b, ODK Error Control 

This task was completed during the reporting period. The 
numerical integration method was rederived, and is attached as 
Appendix A. 

A new error control formula was derived that is based on an 
estimate of the absolute error for the integration step. This 
is done by evaluating the next term (3rd order) in the series 
expansion. The absolute error is converted to a relative error, 
unless the variable is too small, in which case the absolute 
error is used. The error term is used to control step size 
doubling or halving as shown in the attached write-up for 
subroutine INT. 

It was found that accurate integration of the differential 
equations requires consistency in the nozzle cross-sectional 
area and its derivatives. Values of A(x), dA/dx and d 2 A/dx 2 are 
required for the supersonic flow region. The subroutine used 
for interpolation, SPLN, was modified so that y, y * , and y*' are 
found by cubic, parabolic and linear interpolation, 
respectively, such that y* is the integral y ' ' , etc. 

The method was tested using the TDK analysis for the SSME. 
Plots describing the supersonic wall contour and its derivatives 
are shown in Figure 1. Local error vs. x/y* and step size vs. 
x/y* are plotted in Figure 2 and 3- The maximum step size was 
input as .2. The run time for this case was 25 seconds on the 
VAX 780. Previous to implementing the new error control method., 
this case required 127 CPU seconds. The old error control 
forced a fixed step size of .005. Thus, the new method gives a 
factor of 5 improvements on run time for this case. 

The success of this task will prove very useful when task 
3.2a; TDK with ODK Tables, is done since the ODK computer time 
has been reduced substantially. 


nozzle area ratio vs. length ®J nozzle contour 



Figure 1: Plots Describing the Supersonic Nozzle Wall Contour 

(x and y are divided by the throat radius) 
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Local Error vs. Distance 
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Figure 3: Step Size vs. Distance 
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SUBRO UTINE INT 

Provides control for the implicit integration procedure, 
determines the proper set of nonhomogeneous equations to solve, 
and, after each integration step, computes the next integration 
step size according to the following relations: 
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On option, (JF=1) only the fluid dynamic variables are used 
in determining the next integration step size. 

If the step size is halved for the fourth step, the 
integration is restarted using one-half the original step size. 

The correspondence between equation number and physical 
property is: 

Equation Number Property 

1 Velocity of Gas 

2 Density of Gas 

3 Temperature of Gas 

4 NSP+3 Gaseous species mass fraction 

(1...NSP) correspondence to 
( 4 ... NS P + 3) 

When the flow is supersonic, continuity is used to control the 
integration step size to insure that: 


( p v a ) n+ 1 - <p™> n 


(pVA) 


< C0NDEL 


N+1 


where C0NDEL is an input relative criterion with a default value 
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Task 3.2e, Low Temperature Data 

This task has been completed except for adding a default low 
temperature data file to the program. Section 6.1.1 of the TDK 
manual has been rewritten to describe the usage of this data, 
and is attached as Appendix B. Note that only Cp vs T is to be 
input . 

Example, plots obtained from the method are shown in Figures 
4,5, and 6. The properties are for the species H.,. Figure 4 
shows the input table of Cp vs. T. Figures 5 and 6 show the 
integrated values derived from Figure 4 for S and h vs. T. 



Figure 4: C vs T for H, 





Task 3 • 2f , Update JANNAF Data 


The JANNAF thermochemical Data has been updated by obtaining 
a complete new data file from NASA/LRC, and attaching it to the 
program. The letter of transmittal is attached. 


Task 3-2g, Update Kinetic Data 

The reaction rate study has been completed and the results 
were sent to the COR at NAS A/MSFC in September. 


3.3 NOZZLE SHOCK WAVES 

The plan for doing Task 3. 3c, multiple shock waves, has been 
approved by the COR. A small amount of work was done on this 
task during the quarter. 


3-4 VERIFICATION 


No work was done on this task during the reporting period. 


3.5 DOCUMENTATION 


Work was begun on this task near the end of the reporting 
period. 


4.0 WORK PLANNED FOR THE NEXT REPORTING PERIOD 


All remaining work is to be completed during the next 
reporting period. 


5.0 CURRENT PROBLEM AREAS 


The rate of spending on the contract is too low, about 20? 
per quarter rather than 25$ per quarter. It is clear that a 
time extension will be necessary. A revised schedule is being 
pr epar ed . 


6.0 COST STAT US 

As of the end of the reporting period, the cost status of 
this contract is as follows: 

1) Total Cumulative Costs through $63,300 

30 September 1984 

2) Estimated Cost to Complete $55,700 

3) Estimated Percentage of Physical 55$ 

Completion of the Contract 

The costs on the contract to date correspond well with the 
estimated percentage of completion. 
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NUMERICAL INTEGRATION IN ODK 

Given a differential equation of the form d^ = f (x,y), a 

dx 

family of implicit integration schemes is based on the following 
back- di f f er enc i ng differentiation formula. 1 
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For example, using one term in the sum on the left gives 
y „ = y + hy ' which is the usual backward Euler method. 

In general, using the first k terms gives an implicit 
method of order k, with a local error approximately equal to 
the ( k+1 ) term . 

Thus, using the first two terms gives 


or 


7y n.l * 


Vy , - Vy 
J n+ 1 J n 


hy ' 


n+ 1 


| 7y n+1 * 1 7y n 

with a local error approximated by 
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Ref. 1 Dahlquist, G., and Bjorck, A., Nu merical Methods , 
Prentice-Hall, 1974. (see equation 8.3.12). 
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The integration method used in ODK is a variation of the above 
second order method in that y' n is approximated by 
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Solving for k 
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The procedure can be extended easily for a system of n equations 
in n dependent variables. 


Given 


( iy i = f (x,y i , . . ,y n ) , i = 1 , . . , n , 
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equation (1) becomes 
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equation (4) becomes 
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3 k . - 1 k = h ( f + ct . h + Z B . . k. ) 
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and equation (5) becomes 
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The maximum local error is defined by equation (2), 
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APPENDIX B 


6.1.1 THERMODYNAMIC DATA BELOW 300°K. 

If the temperature at any point computed by ODE, ODK, or 
TDK is found to be below the Thermodynamic Data lower 
temperature limit, T^, the polynomial curve fit data (see 
Section 6.1) will be extrapolated to obtain values for the 
thermodynamic data. The extrapolated data can be inaccurate. 
The L0W T CPHS data set can be used to input low temperature 
data for those species for which extrapolation of the JANNAF 
data is not accurate. Only Cp vs. T is input. Enthalpy and 
entropy are obtained by integrating Cp vs. T within the computer 
program . 

The lower temperature limit, T^, in the Thermodynamic Data 
supplied with the program is 300°K. Thermodynamic Data below 
the temperature, T^, may be input by data cards as described in 
Table 6-5. 

An example of this input is given in Table 6-6 which shows a 
card listing extending the Thermodynamic Data for an 0 2 /H 2 
propellant to 100°K. Data in Table 6-6 is taken directly from 
the JANAF tables (Reference 23), except for Argon which is taken 
from NASA SP-3001. 


Ref. 23 Stull, D.R., Prophet, H., et al . , JANAF 
Thermochemical Tables, Second Edition, NSRDS-NBS 
37, National Standard Reference Data Series, 
National Bureau of Standards, June 1971. 


Table 6-5: Input Specifications for Low Temperature 

Thermodynamic Data. 

The first card is a directive card that identifies the start 
of the low temperature thermodynamic data input. It reads as 
follows, columns 1 through 10: 

LOW T CPHS 

The next card contains a species name consisting of 12 
characters, or less, left justified to column 1. An integer is 
placed in column 21 indicating the number of T, Cp values that 
follow. A value of 1,2, or 3 can be used. For example: 

K2 0 2 

Next, the thermodynamic data for the species named above is 
input. There must be one pair of T, Cp values per card. These 
cards are numbered consequently in column 45. They are read as 
2F1 0 . 0 , 20X , 15 . For example, for H^O: 

1 00 . 7.961 1 

200. 7.969 2 

Species name cards and thermodynamic data cards for other 
species, if any, follow. 

A final directive card is used to identify the end of the 
low temperature thermodynamic data. It reads as follows, 
columns 1 through 14: 

END LOW T CPHS 
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TABLE 6 


col.no. 


- 6 . 


LOW TEMPERATURE C° ,H° 
P T T 


S° DATA FOR AN 0 2 /H 2 


PROPELLANT 


I 

I — J.21. ,_45 


LOW T CPHS 


AR 


2 

100,0 

*,9681 


200,0 

*,9681 


H 


2 

100.0 

*.968 


200,0 

*.968 


H2 


2 

100,0 

5.393 


200,0 

6,518 


H 2 o 


2 

100,0 

7.961 


200,0 

7,969 


N2 


2 

100.0 

7 . 07 * 


200,0 

6.989 


0 


2 

100,0 

5.666 


200,0 

5 , * 3 * 


OH 


2 

100,0 

7.567 


200,0 

7,309 


02 

• 

2 

100*0 

6.958 


200*0 

6*961 


END L T 

c*hs 
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1 
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1 

2 

1 

2 

1 
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1 

2 

1 

2 

1 
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